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We investigate a quantum algorithm which simulates efficiently the quantum kicked rotator model, 
a system which displays rich physical properties, and enables to study problems of quantum chaos, 
atomic physics and localization of electrons in solids. The effects of errors in gate operations are 
tested on this algorithm in numerical simulations with up to 20 qubits. In this way various physical 
quantities are investigated. Some of them, such as second moment of probability distribution and 
tunneling transitions through invariant curves are shown to be particularly sensitive to errors. How- 
ever, investigations of the fidelity and Wigner and Husimi distributions show that these physical 
quantities are robust in presence of imperfections. This implies that the algorithm can simulate the 
dynamics of quantum chaos in presence of a moderate amount of noise. 
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I. INTRODUCTION 

It is only recently that it was realized that quantum 
mechanics can be used to process information in funda- 
mentally new ways. In particular, Feynman [1] empha- 
sized that the massive parallelism due to the superpo- 
sition principle may allow to simulate efficiently some 
problems intractable on classical computers, the most 
obvious being many-body quantum systems. Since that 
time, a model of quantum computer has been set up, 
viewed as an ensemble of n qubits, i.e. two-level sys- 
tems, with a Hilbert space of dimension 2" (see reviews 
[2-4]). Computation is performed through unitary trans- 
formations applied to the quantum wave functions of this 
many-body system. In fact, it has been shown that any 
unitary transformation on this 2" dimensional space can 
be written in terms of a set of universal gates, for exam- 
ple one- and two-qubit transformations. Also, important 
quantum algorithms have been developed, such as Shor's 
algorithm for factoring large numbers [5] which is expo- 
nentially faster than any known classical method, and 
Grover's algorithm to search a database [6], where the 
gain is polynomial. 

Motivated by these developments, many experimental 
implementations for actual realization of such a quantum 
computer were proposed (see [4] and references therein). 
Recent results include for example the NMR implemen- 
tation of factorization algorithm with seven qubits made 
from nuclear spins in a molecule [7] , and the simulation of 
the quantum baker map [8]. Thus small quantum com- 
puters with a few qubits are already available experi- 
mentally, and systems of larger size can be envisioned at 
relatively short term. 

Still, algorithms such as the one of Shor require large 
number of qubits and the use of many gates. It is 
therefore important to develop algorithms which need a 
smaller number of qubits and gates and still can yield 
interesting quantities. In particular, algorithms enabling 



to simulate quantum mechanical systems, as originally 
envisioned by Feynman, can be implemented relatively 
easily and solve problems inaccessible to classical com- 
puters with less expanse in number of qubits and gates. 
Several such algorithms have been developed for various 
systems, including many-body Hamiltonians [9] or spin 
lattices [10]. An especially interesting class of systems 
corresponds to chaotic quantum maps. Such systems can 
have a very complex dynamics while their Hamiltonians 
keep a relatively simple form. Algorithms for fast simula- 
tion on a quantum computer were built for the quantum 
baker map [11], the kicked rotator [12], and the sawtooth 
map [13]. We note that recently the quantum baker map 
was implemented on a NMR quantum computer [8] . The 
kicked rotator is an especially rich and generic system, 
which has been a cornerstone for the study of quantum 
chaos [14]. In the classical limit it reduces to the Chirikov 
standard map which has been also extensively studied in 
the field of classical chaos [15]. Implementation of this 
model can be done on a small quantum computer with a 
few tens of qubits, and classical supercomputers will be 
outperformed with a few hundreds of qubits. Still, real 
quantum computers will not be free of imperfections and 
errors, and this will affect the results of the computa- 
tion. It is therefore important to understand the effects 
of different sources of errors on the results of such an 
algorithm. For example, first numerical simulations of 
the quantum computation of this model [16] have shown 
that errors affect in a different way the various physical 
quantities characterizing the model, and that for some of 
them the effect of errors can be exponentially strong. 

In this paper, after presenting in more detail the 
physics of the kicked rotator, we study the effects of er- 
rors on several physical quantities. Wc focus on random 
unitary errors, which may arise when imperfect gates are 
applied, and study first how global quantities such as 
second moment or fidelity are affected by errors. Our re- 
sults confirm and extend those obtained in [16] showing 
a marked contrast in the behavior of these two quanti- 
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tics in presence of errors. We also investigate how well 
the whole wave function is reproduced by an imperfect 
quantum computer. A particularly interesting way to 
display wave functions is to express them through phase 
space distributions, such as Wigner and Husimi func- 
tions. These distributions display the same information 
as the wave functions, but in a form which allows direct 
comparisons between classical and quantum dynamics, a 
property especially interesting to probe the classical limit 
of quantum mechanics. They have been extensively used 
in many fields, and recently a method has been devised 
[17] to measure such distribution for quantum simula- 
tions on quantum computers. The effects of errors on 
Wigner and Husimi functions will be investigated in de- 
tails, showing how imperfections affect the different parts 
of phase space, and discussing how information can be 
retrieved through quantum measurement. A separate 
section is focused on how a localized distribution may 
escape from an island of integrability, showing an espe- 
cially large effect of quantum errors on a quantity which 
is directly relevant to quantum tunneling. 

II. THE KICKED ROTATOR 

The classical kicked rotator is described by the 
Chirikov standard map [14,15]: 

n = n + ksm0; e = e + Tn (1) 

where {n, 9) is the pair of conjugated momentum (action) 
and angle variables, and the bars denote the resulting 
variables after one iteration of the map. It describes a 
free angle rotation and a kick in momentum. This area- 
preserving map has been extensively studied during the 
past decades, and has been applied to problems such as 
particle confinement in magnetic traps, beam dynamics 
in accelerators, comet trajectories and many others [15]. 

The dynamics of this map takes place on a cylinder 
(periodicity in 0), and is controlled by a single parameter 
K = kT. For K = the system is integrable and all tra- 
jectories lie on one-dimensional tori (lines n =constant). 
For K > 0, the system undergoes a transition to chaos, 
which follows the Kolmogorov-Arnold-Moser (KAM) the- 
orem. Periodic orbits corresponding to rational frequen- 
cies are transformed into chains of integrable islands 
mixed with chaotic region. On the contrary, tori with 
irrational frequencies arc deformed but survive, forming 
invariant curves which separate zones in phase space. As 
K is increased, these surviving tori become Cantor sets 
(cantori) and disappear. The most robust torus corre- 
sponds to the golden number (l-|--\/5)/2, and disappears 
for K = Kg7^ 0.9716.... Thus, for K > Kg global chaos 
sets in, with appearance of an extended chaotic region 
in phase space, and with dynamics characterized by a 
positive Kolmogorov-Sinai entropy h « \q.{K/2) > (for 



K > 6). In this regime, a typical trajectory shows dif- 
fusive growth of momentum which statistically can be 
described by the Fokker-Planck equation, with diffusion 
rate D = -n? /t « where t is measured in number of 
iterations (kicks) [14,15]. For lower values of iiT, the phase 
space displays a complex hierarchical structure with in- 
tegrable islands surrounded by chaotic zones at smaller 
and smaller scales. 

The map (1) is periodic in n with period 2-k/T, so 
the phase space structures repeat themselves on each cell 
of size 27r/T. Such a cell is shown on Fig.l for K = 
Kg, displaying the complex hierarchical structures which 
appear in the phase space. 




FIG. 1. Plot of the classical phase space at K = 
Kg = 0.9716... {t = 10* iterations of (1) for 200 points). 

The quantization of (1) yields a Hamiltonian which 

after integration over one period gives a unitary evolution 
operator acting on the wave function -0 

^ = f7V' = e-*'='=°"^e-*^"'/V, (2) 

where h = -id/dO, h = 1, and ■0(6' + 27r) = ■(/>((?). The 
quantum dynamics depends on two parameters k and T 
(instead of the single parameter K = kT for the classical 
one). The classical limit corresponds to fc oo, T — > 
while keeping K — kT =constant [14,18,19]. In a sense, 
T plays the role of an effective %. 

Depending on the values of theses parameters, the sys- 
tem follows different regimes, from regularity to quantum 
chaos. Due to this variety of behaviors, the quantum 
kicked rotator has been intensively studied (see [14,18,19] 
and Refs. therein). Indeed, most of the phenomena char- 
acteristic of quantum chaos are present, such as quantum 
ergodicity. Random Matrix Theory statistics, chaos as- 
sisted tunneling, and others. In particular, for K > Kg, 
the phenomenon of dynamical localization appears. Al- 
though in this regime a typical classical trajectory dif- 
fuses in momentum, the eigenstates Xm{n) of the opera- 
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tor U in momentum space arc exponentially localized for 
typical values of k and T. Their envelopes obey the law 
Xmin) ~ exp(— |n — m\/l)/Vl where m marks the cen- 
ter of the cigcnstatc and / is the localization length. For 
k ^ K 1 this length is determined by the classical dif- 
fusion rate I = D/2 w A;^/4 [18]. This phenomenon has 
close relationship with the Anderson localization of elec- 
trons in disordered solids [20], and investigation of the 
kicked rotator gives information on this important solid- 
state problem still under intensive investigation nowa- 
days. The quantum kicked rotator describes also the 
properties of microwave ionization of Rydberg atoms [21]. 
It has been realized experimentally with cold atoms, and 
the effects of dynamical localization, external noise and 
decoherence have been studied experimentally [22]. 

For numerical studies of the quantum evolution (2) it 
is convenient to choose the case of quantum resonance 
with T/(47r) = M/N where M,N are integers [19]. In 
this way the quantum dynamics takes place on a torus 
with N levels. For I ^ N the eigenstates of evolution 
operator become ergodic and the level spacing statistics 
is described by random matrix theory [19]. 

The algorithm for the quantum simulation of the 
kicked rotator was presented in [12]. The evolution (2) 
consists of the product of two unitary operators, which 
are diagonal in the angle and momentum bases respec- 
tively. The most efficient classical algorithm available 
consists in changing back and forth between the angle 
and momentum representation by Fast Fourier Trans- 
forms (FFT) . The operator which is diagonal in the basis 
is then implemented by direct multiplication of the co- 
efficients of the wave function. In total, one iteration of 
(2) on a Hilbert space of dimension N = 2"'' requires 
O(A^logA^) classical operations, the limiting steps being 
the FFT. The quantum algorithm follows the classical 
one, and speeds up all parts of it to obtain exponen- 
tial increase of compiitation rate. First an initial dis- 
tribution is built, in a polynomial number of operations 
(in Hq). Various initial wave functions can be built in 
such a way. In the following, we will use as initial state 
l^fo) a- wave function localized at a precise value of mo- 
mentum n, which can be built in Uq single-qubit rota- 
tions starting from the ground state. The general state 
of the system can be written as ^^n^Qdnln), where a„ 
are the amplitudes of the wave function on the \n) ba- 
sis state. Then the first unitary operator is applied. In 
the n representation it is diagonal and can be written 
exp(— iTn^/2). This operator can be implemented effi- 
ciently by using the binary decomposition of n: if n = 
Ejlo'"j2^ then n2 = Ej„j,aj,aj,2^'+^\ Therefore 
exp(-iTn2/2) = Uj^j^exp{-iTaj^aj^2^^+^''-'^) with 
'^ji 2 = or 1 . Thus one needs to implement the two- 
qubit gate applied to each qubit pair (ji, j2) which keeps 
the states |00), |01), |10) unchanged while |11) is trans- 
formed to exp(-ir2-'i+-?2-i)|ll). O(n^) applications of 



this gate arc sufficient to simulate cxp(— irn^/2). 

Then a quantum Fourier transform (QFT) (see e.g. [2]) 
is performed to shift from n to ^ representation, yield- 
ing YliLa^ ■ This transformation needs only 0{n'^) 
one and two-qubit gates, and yields the wave function 
in 6 representation. In this representation, the second 
operator cxp(— i/ccos^) is diagonal. Direct (sequential) 
multiplication by exp(— ifccos^j) for each 0i will require 
exponentially many operations, so a parallel way to ap- 
ply this operator has to be devised. In [12], it was pro- 
posed to use supplementary registers on which the values 
of cos(^i) will be computed in parallel. The procedure 
transforms I^j^q^ &i|6'i)|0) into Y^fj^^ h\9i)\ cosOi) with 
cos{0i) computed up to a fixed precision using a recursive 
method based on Moivre's formula [12]. This is actually 
the slowest step of the algorithm, requiring 0{n^) ele- 
mentary operations. From the state '^^Jq^ bi\Si)\ cosOi), 
it is easy by using nq one-qubit operations to build 
the state ^^^^^ 6i exp(— ifccos^i)]^i)| cos0,). Then the 
cosines in the last register are reversibly erased by 
running backward the sequence of gates that con- 
structed them, and one ends up with the state 
Y^f^^ bi cxp(— ik cos 9i)\9i)\0), which is the result of the 
action of the unitary operator exp{—ik cosO). Another 
QFT (requiring O(n^) operations) takes the wave func- 
tion back to the n representation. 

Increasing Uq, which exponentially increases the di- 
mension of the Hilbert space available, enables to probe 
various physical limits in the system. If K = kT is kept 
constant, the classical mechanics remains the same. If 
T is kept constant, the effective % is fixed, and increas- 
ing Uq will increase exponentially the size of the phase 
space of the system (number of cells). In contrast, if 
T = 2'k/N, with AT = 2"«, the size of the phase space 
remains the same, all A'^ momentum states correspond- 
ing to the same cell of size 2tt/T. In this case, increasing 
Uq increases the number of quantum levels correspond- 
ing to the same classical structure, and is equivalent to 
decreasing h toward the classical limit. 

The whole quantum algorithm described above re- 
quires 0(?T.q) gate operations to perform one iteration of 
the quantum map (2), exponentially less than the clas- 
sical algorithm. Still, a physical quantum computer will 
not be an ideal perfect machine, and there will be im- 
perfections, which may hamper the computation. In the 
following sections, we will investigate the effects of noise 
and imperfections on the physical quantities that are sim- 
ulated, and estimate the accuracy of the quantum com- 
putation of the quantum map (2). The numerical sim- 
ulation of many qubits is very resource-consuming on a 
classical computer. Due to that we took in all numerical 
computations the action of exp(— ifccos^) as exact, and 
performed by direct multiplication in 6 representation, 
all other operations being made with errors. We think 
that this approximation does not alter the qualitative 
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features of the results, although the number of quantum 
gates is reduced from O(n^) to O(n^). Also in this ap- 
proximation all supplementary registers required for the 
computation of | cos^^) are eliminated and the quantum 
evolution on iV = 2^1 levels is performed with only Uq 
qubits. 

III. GLOBAL QUANTITIES 

We first study the effects of imperfections and errors 
for the global quantities of the system. 
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FIG. 2. Dependence of seeond moment {v?) — {{n — no)^) 
of the probability distribution on time t for T = 0.5 
and K = 15. Data are shown from top to bottom for 
Uq = 16,15,14,13 and e = 10^'' (four curves). The low- 
est fifth full curve is for e = 0, n, = 14. Initial state is 
|*o> = |no>, with no = N/2. 

To model these imperfections, we introduce a random 

unitary error during the operation of elementary quan- 
tum gates. These errors are present for each gate per- 
forming the quantum Fourier transform and the action of 
the unitary operator e^*"^" Two elementary gates arc 
used: single-qubit Hadamard gates H = diag(l, 1,1,-1) 
and the two-qubit gate B = diag(l, 1, l,exp(ia)) where 
a is a phase. The transformation H can be written H = 
Uqct where uo = (l/v^, 0, l/-\/2) and a = {ax, cry, Uz)- It 
is replaced by an imperfect gate H' = ua, where u is a 
unit vector with a random angle (3 from uq. In a, similar 
way, each B is replaced by B' = diag(l, 1, 1, exp(ia-|-i7)) 
where 7 is again a random angle. At a given strength 
£ > of noise, each gate is implemented with a /3 or 
7 randomly selected from a uniform distribution such 
that < 7re or I7I < tts [23]. As explained in Section 
II, we made the approximation of taking the action of 
exp(— ifccos^) as exact, all other operations being made 
with errors. The use of Uq qubits gives a Hilbert space 
for wave functions of the kicked rotator on N levels, with 
N = 2"" , i. e. values of momentum range from n = 1 to 
n = N. In all numerical computations, the initial state 



l^'o) was chosen as localized on a precise value of the mo- 
mentum no, i. e. I^'o) = |n.o), with no = 1 (lowest value 
of momentum) or no = A''/2. The rotation is computed 
as exp(-iT(n - n)2/2) with n = N/2. 

Depending on the choice of parameters in (2), increas- 
ing the number of qubits rig will increase the number 
of values of momentum in each phase space cell of size 
An = 2tt/T, or increase the number of cells, or both. In 
[16], it was shown that if T is constant while n, increase, 
errors in the QFT may lead to an exponential growth of 
errors with Uq for the second moment (n^) of the prob- 
ability distribution. In this case, the size of phase space 
grows exponentially with n, but K and the effective h 
are kept fixed. Due to quantum localization, exact wave 
functions cannot spread beyond a region of size given 
by the localization length, which remains fixed when Uq 
increases. Therefore for all values of Uq, the second mo- 
ment of a distribution initially located at no = A^/2 will 
saturate with time at a value independent of Uq (full line 
in Fig. 2) if (2) is exactly simulated. On the contrary, er- 
rors in the QFT lead to small transfer of probabilities to 
regions of phase space which are exponentially far away 
from where the exact wave function is localized. This 
induces the exponential increase of the second moment 
with Uq. We confirm here this effect in Fig. 2 for different 
parameters with the more complete set of errors used in 
this paper, and with simulations up to larger number of 
qubits. 
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FIG. 3. Dependence of rescaled time scale tq on number 
of qubits Hq for 10"® < e < 0.03, T = 0.5, K = 5 (x) and 
K = 15 (o). Initial state is |5'o) — \no) with no = N/2. Data 
are averaged over 10 to 1000 reahzations of noise. Full and 
dashed lines correspond to the theoretical formula (3) with 
Cq = 0.23. Logarithm is decimal. 

To be more quantitative. Fig. 3 shows the time scale tq 
on which the presence of errors leads to a doubling of the 
value of the second moment (n^) as a function of Uq and 
error strength e. In [16] the formula: 

tq^Cqk*/{e^nq2^^^) (3) 
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was proposed and checked numerically with up to 13 
qubits. It stems from the fact that each imperfect gate 
operation transfers on average a probability of equally 
divided among Uq spurious peaks located at integer pow- 
ers of 2. Thus due to imperfections (n^) ^ nqe'^2'^"-''t 
(each time step involves ^ gate operations), whereas 
for the exact wave function (n^) w w 4P w 
Both expressions become comparable at the time tg given 
by (3). Fig. 3 confirms this formula by extensive numeri- 
cal computations, with up to 20 qubits, and for two dif- 
ferent values of K. This enables to get the numerical 
constant Cq « 0.23. 




FIG. 4. Evolution of fidelity / with time t. Full curves are 
for K = 1.3, T = 2n/N {N = 2""), Uq = 14, and from top 
to bottom e = 3 X 10"^, e = 10"^, e = 0.03. Initial state 
is l^'o) = |no) with no = 1. Dashed curves are for T = 0.5, 
K = 5 and Uq = 14, and from top to bottom e = 3 x 10~^, 
£ = 10~^. Initial state is |>]^o) = |^^o) with no = N/2. Loga- 
rithm is decimal. 

Although the time scale tq drops exponentially with Uq, 
there are other observables which show only polynomial 
sensitivity to errors. A standard quantity used to charac- 
terize the global influence of errors is the fidelity defined 
by the projection of the wave function with errors i/'e (t) 
on the perfect one ^Q{t): f{t) = \{il)^{t)\i>o{t))\'^ . The 
dependence of this fidelity on time in presence of errors 
is shown in Fig. 4, showing that it decreases slowly with 
t and amplitude of noise e. One can define a time scale 
tf such that f{tf ) = 0.5. Fig. 5 presents the variation of 
tf with system parameters in two different regimes. It 
shows that the relation: 

tf « Cf/ie^v:^,) (4) 

holds with the numerical constant Cy ~ 0.35. Fig. 4 and 
Fig. 5 are consistent with a fidelity decay f{t) ~ exp(— Fi) 
where F ~ e^n^. 

The relation (4) can be understood from the following 
physical considerations. Each imperfect unitary gate is 



rotated by a random angle of order e from the exact one. 
Therefore a probability of order is transferred from 
the exact state at each gate operation. Each time step 
of the map (2) takes 0(?^q) operations, in the approxi- 
mation which we have taken where the building of the 
cosines is supposed exact. This implies that tf, which is 
in units of time steps of (2), should vary as 1/ (e^n^). We 
expect that if the full algorithm was implemented, with 
the cosines computed following the procedure explained 
in Section II, a time step of (2) should take O(n^) oper- 
ations, and accordingly tf should vary as l/(£^nq). 
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FIG. 5. Dependence of time scale tf on system parameters 
for n, = 4 (o), 6 (□), 8 (o),10 (A), 12 (<), 14 (v), 16 (>), 18 
{+). Here K = 1.3, T = 2tt/N [N = 2"") (open symbols) or 
K = 5, T = 0.5 (full symbols). Dashed line is the theoretical 
formula (4) with C/ — 0.35. Initial state is I'I'o) = |no), with 
no = 1 {K = 1.3) or no = N/2 {K = 5). Data are averaged 
over 10 to 100 realizations of noise. Data for T = 0.5 and 
K = 15 are nearly indistinguishable from T = 0.5, K = 5 
(not shown). Logarithms are decimal. 

The data shown in this section exemplify the sharp 
contrast in the behavior of the different observables in 
presence of errors. The fidelity shows only a polyno- 
mial decrease with respect to both e and Uq, whereas 
the second moment of the wave function grows exponen- 
tially with Uq, but polynomially with e. The resolution of 
this apparent paradox is related to the fact that the sec- 
ond moment is sensitive to the size of the Hilbert space, 
which grows exponentially with Uq . Small spurious peaks 
due to imperfections do not spoil the fidelity, but modiiy 
strongly the variance (n^) if they appear very far away 
from the exact location of the wave function [24]. 
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IV. WIGNER AND HUSIMI DISTRIBUTIONS 

In the previous section we focused mainly on the case 
where T (effective fi) is fixed but the phase space size 
grows exponentially with Uq. In contrast, at T = 2-k /N 
and N — 2"' the system size in classical momentum 
(number of 2tt/T cells in n) remains fixed when rig in- 
creases. In this way the effective Ti drops exponentially 
with Uq and going to larger number of qubits means ap- 
proaching the classical limit (exponentially fast). Smaller 
and smaller details of the classical structure will be visi- 
ble in the quantum wave functions. In this regime, data 
presented in Fig. 4 and Fig. 5 have already shown that the 
fidelity follows the law (4) as in the case T =constant. 
However, the fidelity characterizes in one number the ac- 
curacy of the whole wave function, and does not tell how 
well the local properties are reproduced. To study the 
local properties of wave functions, one can express it in 
6 OT n representation. However, a very useful represen- 
tation corresponds to phase space distributions, such as 
Wigner or Husimi distributions. They are especially used 
in the field of quantum chaos, since such representations 
permit a direct comparison with classical Hamiltonian 
mechanics which takes place in phase space. They also 
enable to probe the classical/quantum border when h is 
decreased compare to other parameters of the system. 
Plotting such quantities in presence of errors allows to 
probe how local properties of the wave functions are sen- 
sitive to imperfections in the quantum algorithm. 

An additional motivation to study such phase space 
representations stems from the fact that recently an al- 
gorithm was proposed [17] which enables to compute the 
Wigner function on a chosen point in phase space by the 
use of an ancilla qubit. 

For a continuous system with two conjugate variables 
p and q the Wigner transform [25] of a wave function tp 
is defined by: 



W{p,q) 



-^{q+\yiK<l-\)dq' (5) 



In a discrete system with iV-dimensional Hilbert space, 
one is led to define the Wigner function on a lattice of 
2A^ X 2A^ points (see e. g. [26]). In the case of the kicked 
rotator, the formula becomes: 



N-l 
m=0 

NB 



'-n{m,-e/2) 



2N 



-iP{e-m)*i;{m), (6) 



with = ■j;^. The Wigner function is always real, but 
contrary to classical Liouville phase space distributions 
it can take negative values. It verifies J^i^i^i^''^) ^ 
|V(n)|2 and J:^Wi0,n,) = ■ 




FIG. 6. (color) Plot of Husimi (left) and Wigner (right) 
distributions at t = 10^ for K ^ 1.3 > Kg, T = 2tt/N, 
N = 2"' and Uq — 7. Initial state is 1^*0) = I'^o), with 
no — 1. Top; e = 0; middle: e = 0.002; bottom: 
£ — 0.004. Left: color/grayness represents intensity level from 
blue/white (minimal) to red/black (maximal). Right: gray- 
ness represents amplitude of the Wigner function, from white 
(minimal negative value) to black (maximal positive value). 

The Wigner transform has the drawback of being 
negative or positive. Nevertheless, coarse-graining this 
function over cells of size h gives non negative values. 
Such a procedure gives the Husimi distribution (see e.g. 
[27]) which corresponds to a Gaussian smoothing of the 
Wigner function. In the case of the kicked rotator, the 
Husimi distribution can be computed through: 



h{9,n) 



n+N/2 

E 

m=n-N/2 



N 



(7) 



where the Gaussian for simplicity is truncated for val- 
ues larger than N/2, and ^p{m) is the wave function in 
momentum representation. The Husimi distribution is 
always non negative, and allows a direct comparison be- 
tween classical Liouville density distributions and quan- 
tum wave functions. 

Wigner and Husimi distributions of wave functions of 
the quantum kicked rotator simulated on a quantum com- 
puter are shown on Fig. 6 for different level of errors. Both 
functions have similar patterns, although as expected the 
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Wigner function displays interference structures absent 
in the Husimi distribution. In the regime of parame- 
ters studied, classical invariant curves are still present in 
phase space and prevent the exact wave function to enter 
the large elliptical island in the middle. In the presence of 
moderate level of noise, main structures are still present 
and distinguishable. 




FIG. 7. (color) First three rows: Husimi distribution at 
t = 10^ for if = 1.3 > Kg and T = 2tt/N, N = 2""; from 
top to bottom: Uq = 9, Uq — 12, Uq — 14; quantum noise 
£ = (left), e = 0.002 (center), e = 0.004 (right). Bot- 
tom row: classical phase space distribution at t = 10^ with 
classical noise £ = (left), e — 0.002 (center), e — 0.004 
(right). For clarity, the distributions are averaged over 10 
iterations around t — 10^. Initial quantum/classical state 
is no = 1. Color/gray ness represents intensity level from 
blue/white (minimal) to red/black (maximal). 

Fig. 7 confirms this result, showing the Husimi dis- 
tribution for larger number of qubits, together with the 
classical phase space distribution. The Husimi distribu- 
tions in phase space show features mimicking the classical 
phase space distributions, in accordance with the corre- 
spondence principle. Fig. 7 (left) shows that when Ug is 
changed, finer and finer details of the classical structures 
are visible in the exact quantum wave function, in accor- 
dance with the fact that increasing Uq amounts to reduce 
h, and approach the classical limit. The same figure shows 



that the wave function is spread over a larger domain of 
phase space as Uq increases. This can be explained by the 
following effect. In this mixed regime between integra- 
bility and strong chaos at i^T = 1.3 the invariant classical 
curves which prevent any transport arc no longer present 
since the last one is destroyed at K = Kg = 0.97... But 
cantori are present, which are remnants of the disap- 
peared invariant curves. They have a fractal structure, 
and a wave packet can cross them only if the holes are 
large enough. These holes scale as (K — Kg)^ and become 
comparable with the minimal area scale of the Husimi 
distribution determined by the effective h given by T. 
Hence, for K — Kg ^ 1, the wave function is prevented 
to cross the cantorus for {K — KgY < T. Due to that 
quantum interference prevents the transport via cantori 
[14,18,28]. 

The quantum Husimi distributions shown in Fig. 6 and 
Fig. 7 display structures of increasing complexity with 
larger nq. Still, with moderate level of noise, the quan- 
tum computer is able to reproduce the exact distribu- 
tions with reasonable accuracy. For larger errors in gate 
operations, significant probability is present at wrong 
phase space locations, and phase space structures be- 
come blurred. The comparison with the effect of classical 
noise visible in Fig. 7 shows that in this case the quantum 
errors enable the wave function to enter classically for- 
bidden zones much faster, a fact which will be analyzed 
in more details in the next section. 

0.8 I ' < ■ < ■ < ■ < ■ 1 
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FIG. 8. Relative error on the 

Wigner function SWe — {\W — We\) / (\W\) as a function of 
time ior K = Kg, T = 2n/N, N = 2"' and Uq = 10. Ini- 
tial state is I'I'o) — \no), with no = N/2. From bottom to 
top quantum noise is £ = 10~^, £ = 10"'^'^, £ — 10"'^. The 
Wigner function is averaged over 2A'' values in the chaotic 
zone. Data are averaged over 10 realizations of noise. 

It is interesting to evaluate the effects of noise and 
imperfections not only on the broad features of the full 
Wigner function, but also on individual values. In Fig. 8 
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and Fig. 9, the behavior of individual values of the Wigncr 
function in presence of noise in the gates is investigated. 
Fig. 8 shows that the relative error (i.e. the error (IW — 
W^\) divided by the average individual value of the exact 
Wigner function increases slowly with the growth 

of t and e even in the chaotic zone. Similar results can 
be observed in the intcgrablc zone and in the localized 
regime (data not shown). In a more quantitative way, 
Fig. 9 shows the behavior of the time scale tw when the 
error on the Wigncr function become comparable to its 
mean value in the regime chosen ((|M^(fH') — l^e(iw)|) = 
(|T4^|)/2). In all three cases considered, one obtains 

tw « Cw/iW^e^), (8) 

with a = lora = 1.5. Thus individual values of 
the Wigner function in the kicked rotator model are ro- 
bust quantities with respect to noise, even in the chaotic 
regime. These results are interesting also in view of the 
recent discussion on the effects of decoherence on Wigner 
functions [29] . Our results clearly show that in the frame- 
work of quantum computation, the errors on the Wigner 
function are polynomial and not exponential. 
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log(nq 8) 

FIG. 9. Dependence of time scale tw on system parameters 
for 5 < < 11. Here K = Kg,T = 27r/7V (JV = 2"'). The 
Wigner function is averaged over 2N values in the chaotic 
zone (O) or in the intcgrablc zone (A). Straight lines arc 
the theoretical formula (8) with a = 1.5 and Cw = 0.02 (full 
line) or Cw = 0.03 (dashed line). Initial state is Y^o) = |no), 
with no = N/2. Data are averaged over 10 to 1000 realiza- 
tions of noise. Inset: Dependence of time scale tw on system 
parameters for 5 < riq < 14. Hero T = 0.5 and K = 5. The 
Wigner function is averaged over 2A'" values in the localized 
zone (n). Full line is the theoretical formula (8) with a = 1 
and Cw = 0.012. Initial state is [*()) = [no), with no = N/2. 
Data are averaged over 10 to 1000 realizations of noise. Log- 
arithms are decimal. 

As noted previously, a recent algorithm [17] enables to 



measure the value of the Wigner function of a system of 
density matrix p on a selected point in phase space, with 
the help of an ancilla qubit a. First H (Hadamard gate) 
is applied on a, followed by a controlled-?/ operation (U 
is applied to the system to be measured depending on 
the state of a) and again H is applied on a. Then the 
expectation value of a is (cr^) = Re[Tr{U p)]. The use 
of a particular operator U, which can be implemented 
efficiently [17], enables to get W{p,q) = {a^)/2N (where 
A/' = 2"<'). 
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FIG. 10. Dependence of inverse participation ratio ^ of 
the Wigner function on number of qubits riq at t = 10^ for 
T = 2ti/N, N = 2"9 and K = 0.5 (full curve), K = 0.9 
(dashed curve), K = 1.3 (long-dashed curve), K = 2.0 
(dot-dashed curve). Initial state is |*o) = \no), with no = 1, 
and e = 0. 
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FIG. 11. Dependence of inverse participation ratio ^ of 
the Wigner function on number of qubits Uq at t = 10^ for 
T = 27r/iV, N = 2"" and K = 2 (full line), and T = 0.5 
and K = 5 (dashed line). Dotted lines show ^ oc iV^ and 
^ oc iV. Initial state is j^to) = jno), with no = 1, and £ = 0. 
Logarithm is decimal. 
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However, we should note that even if U can be im- 
plemented efficiently, (cr^) can be evaluated only by it- 
erating the procedure enough times to get a good esti- 
mate. Therefore the amplitude of the signal is crucial to 
make the whole process efficient. Thus it is interesting 
to study the amplitude of peaks in the Wigner function, 
in order to know if strong peaks are present which can 
be detected reliably through this method. This can be 
investigated through a quantity which we call inverse par- 
ticipation ratio of the Wigner function, in analogy with 
the inverse participation ratio for wave functions used in 
quantum chaos and systems with Anderson localization 
[30]. For a wave function with N projections ipi on some 
basis, the inverse participation ratio X) IV'il^) 
measures the number of significant components in this 
basis. For the Wigner function, one has the additional 
sum rules ^ VFj = 1 and J2 — W- '^'^ define an 
inverse participation ratio for the Wigner function, we 
therefore use the formula ^ = 1/(7V^ ^ Wj^)- K ^ peaks 
of approximately equal weights 1/7V are present, then 
^ = N, whereas iV^ components of equal weights (in 
absolute value) give ^ = N"^. This quantity ^ 

therefore permits to estimate the number of main com- 
ponents of the Wigner function. Fig. 10 and Fig. 11 show 
the scaling of this quantity with Uq for different values of 
parameters. In all the cases where T = 2'it/N {N = 2"') 
the ratio ^/N'^ reaches a saturation value. This implies 
that asymptotically (ctz) = NW{p,q) ~ 1/VN, a value 
which requires A'^ iterations followed by measurements to 
be reliably estimated. In this case, the asymptotic gain in 
number of operations compared to the classical algorithm 
is only 0(log(A^)), although the resources needed are ex- 
ponentially smaller (uq qubits instead of 2"« classical reg- 
isters). This should be contrasted with the case where 
the number of cells increases (T constant) where Fig. 11 
shows that ^ ~ A". This gives (ct^) = NW{p, q) ~ 1, 
which means that in this regime with localization, any 
of the ^ N components of the Wigner function which 
arc important can be estimated reliably and efficiently 
through this method (provided one knows beforehand the 
approximate position of the localized state). The results 
presented in Fig. 8 and Fig. 9 show that despite the dif- 
ferent scaling laws of the relative errors grow only 
polynomially in all cases considered, thus enabling such 
measurements of individual values of W to be reliable for 
moderate amounts of noise. 



V. QUANTUM TUNNELING THROUGH 
INVARIANT CURVES 

In the previous section, it was shown that the classi- 
cal and quantum errors affect the dynamics in a rather 
different way. This difference is particularly striking in 
the regime where classical invariant curves are present 
(integrable or mixed systems, which correspond to mod- 



erate values oi K here, as in Fig. 7). Such invariant curves 
cannot be crossed classically, and only quantum tunnel- 
ing can transfer probability inside integrable islands from 
chaotic regions. However, whereas small classical errors 
enable to cross only neighboring invariant curves, small 
quantum errors may lead to long-distance "jumps" of 
probability deep into integrable island (see Fig. 7, last 
column) . 
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FIG. 12. Dependence of the probability I of the Husimi 
distribution inside the circle (see text and inset) on time t for 
£ = 10-3 and ng = Uat K = 1.3 and T = 2n/N {N = 2"'). 
Initial state is |^'o) = |?^o), with no — 1. Inset: position of 100 
points initially at no = 1 after 10'' iterations of the classical 
map (1), and location of the circular domain D (see text). 
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FIG. 13. Dependence of the probability I of the Husimi 
distribution inside the circle (see text and Fig. 12) on for 
K = 1.3 and T = 2n/N, N = 2"\ s = 3 x 10"^ (solid 
curve) and e = (dashed curve). Data arc averaged over 100 
iterations around t = 10^. Initial state is |4'o) = |no) with 
no = 1. Logarithm is decimal. 
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To study the effects of errors on quantum tunneling we 
show in Fig. 12 the dependence of probabihty of Husimi 
distribution h{6, n) inside the classically forbidden region 
on time t. The quantity I(t) = Jj^ h(9,n)d0dn, where D 
is the domain enclosed by the circle in Fig. 10 (inset), 
shows a linear growth with t. This can be understood 
by a physical argument similar to the one justifying (4). 
Indeed, imperfect gates transfer on average a probabil- 
ity of order from the exact wave function to wrong 
phase space positions. However, not all gates will trans- 
fer probability inside D but only a subset of them. This 
predicts that I{t) rigS^t. Data from Fig. 12 and Fig. 13 
and additional data (not shown) confirm this prediction, 
with a ~ 1.3. 

To exemplify the effect of quantum errors. Fig. 13 shows 

I{t) at fixed time i as a function of number of qubits 
Uq for zero and nonzero noise in the gates. In the case 
of zero noise, there is an exponential decrease with rig. 
Indeed, the only process which allows to enter the island 
for the wave packet is quantum tunneling. In general, 
the probability of such a transition scales like exp{—S/Ji) 
where S' is a classical action. Increase of Uq amounts to 
decrease the effective h and leads to the exponential drop 
of I obtained numerically at £ = 0. In sharp contrast, 
the presence of imperfections in the gates {e > 0) leads 
to direct jumps inside the island that gives an increase 
of I with riq according to the estimate of the previous 
paragraph. Thus for this specific process, the effect of 
noise in the gates results in a qualitative change of the 
dependence of tunneling probability 7 on n,. 



VI. CONCLUSION 

The results presented in this paper show that it is pos- 
sible to simulate efficiently the quantum kicked rotator on 
a quantum computer. For the quantum algorithm simu- 
lating the dynamics of kicked rotator we investigated the 
effects of gate errors and showed that certain quantities 
like fidelity, Wigncr and Husimi distrilnitions arc suffi- 
ciently robust against noise in the gates. Thus for small 
amplitude of noise these quantities can be computed re- 
liably without application of quantum error corrections. 
At the same time we found that there exist other char- 
acteristics, e.g. variance of probability distribution and 
tunneling pro})ability inside stability islands, which are 
very sensitive to errors in quantum gates. In addition, 
the study of the Wigner function shows that individual 
values of this function are robust with respect to quantum 
errors and can be reliably estimated. However, the com- 
putation of the Wigner function at specific points meets 
certain readout problems in deep quasiclassical regime 
where generally a large number of measurements is re- 
quired. 

On the basis of obtained results we believe that the 
quantum algorithms simulating quantum chaotic maps 



will provide important grounds for testing the accuracy 
of the next generation of experimental implementations 
of quantum computers. 
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